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Abstract 

Using fully minimized fundamental measure functionals, we investigate free energies, vacancy 



^4— I ■ concentrations and density distributions for bcc, fee and hep hard-sphere crystals. Results are 

O ; 

. . complemented by an approach due to Stillinger which is based on expanding the crystal partition 



function in terms of the number n of free particles while the remaining particles are frozen at 



B 

T^ • their ideal lattice positions. The free energies of fcc/hcp and one branch of bcc agree well with 

Ch ■ 

Q . Stillinger's approach truncated at n = 2. A second branch of bcc solutions features rather spread- 

o ■ 

out density distributions around lattice sites and large equilibrium vacancy concentrations and is 
presumably linked to the shear instability of the bcc phase. Within fundamental measure theory 
and the Stillinger approach (n = 2), hep is more stable than fee by a free energy per particle of 

o. 

V^ , about 0.001 k^T. In previous simulation work, the reverse situation has been found which can be 



rationalized in terms of effects due to a correlated motion of at least 5 particles in the Stillinger 
picture. 



I. INTRODUCTION 

The crystal lattices of monatomic substances are very often of face-centered cubic (fee), 
hexagonally close-packed (hep) or body-centered cubic (bcc) type. Still, it is a formidable 
problem in statistical mechanics and quantum chemistry to predict the stable crystal struc- 
ture and its free energy for a given substance. Approximating the particle interactions in 
this substance by classical two-body potentials makes the problem amenable to a treatment 
using methods of classical statistical mechanics, most notably Monte Carlo (MC) simulations 
and (classical) density functional theory (DFT). While the approximation using two-body 
potentials may not be very accurate for truly atomic substances, the advance in colloid 
synthesis allows to realize systems with simple two-body potentials to a good degree of ap- 
proximation, thus colloid suspensions are a perfect model system for investigating freezing 
in classical statistical mechanics. 

For isotropic two-body potentials u{r) (r is the center distance between two particles) 
a substantial amount of knowledge has been gathered. For potentials with a repulsive core 
the steepness of the core mainly determines the stability of fee over bcc, with fee being more 
stable for steeper cores. This has been investigated for power-law potentials u oc (1/r)"' 

n Tin 

[ij and screened exponentials u oc exp{—Kr)/r [2, l3| where the parameters n, k determine 
the steepness of the potential. In the hard-sphere limit (n, k — > oo), fee appears to be the 
stable, equilibrium structure and a possible bcc structure is unstable against small shear 
^ which is reflected in squared phonon frequencies a;^(k) being negative for certain wave 
vectors k. 

For hard spheres, it is a much more delicate issue whether fee is more stable than other 
close-packing structures, most notably hep. Early theoretical work by Stillinger et al. an- 
alyzed the free energy of hard disks and fee and hep hard sphere crystals in terms of an 
expansion in the number n of contiguous particles (free to move) in an otherwise frozen 
matrix of particles at their ideal lattice positions 5|-[7| (see below). This expansion could 
be done analytically only for densities in the vicinity of close-packing and, for n = 2 and 
n = 3 (by quite a tour de force), resulted in hep being more stable than fee by a free energy 
difference per particle AF/N ~ 10"'^ k-sT. However, the individual terms contributing in 
this series are much larger than this value of AF/N. An extension of this method [8| (still 
only near close-packing) to n = 5 shows the reverse situation: fee is more stable than hep 



and AF/N ~ — 10~^ k^T, but the last term in the series is still larger in magnitude than 
AF/N (about 6 times for fee and 3 times for hep). Simulation work confirms the stability 
of fee over hep also for smaller densities (around coexistence). Using a single-occupancy cell 



I9 



(SOC) method, Ref. J9| estimates AF/N = -(5± 1) ■ 10^^ ksT at a density of poa^ = 1.041 
(approximately at coexistence, a is the hard sphere diameter). In this method, particles 
are constrained to their Wigner-Seitz cells and the free energy difference is found by in- 
tegrating the equation of state. The limitations of this method could be overcome by the 



lows to compute directly the 
10] . At pocr^ = 1.10 the result 



powerful Monte-Carlo (MC) lattice switch method which a 
free energy difference between two different lattice structures 
is AF/N = -(0.86 ± 0.03) ■ 10^^ AjbT. Thus the result of the high-density Stillinger series 
for n = 5 for the stability of fee over hep and the magnitude of the free energy difference 
is consistent with the MC simulation result at a considerably smaller density. One may 
tentatively conclude that for all densities the stability of fee in the hard sphere system is a 
subtle result of the correlated movement of five and more particles and the effect in the free 
energy is very small. 

In view of this evidence it appears to be very hard to contribute to the theoretical un- 
derstanding of the stability of fee over hep beyond the Stillinger arguments. In this respect, 
density functional theory (DFT) seems to be the only promising candidate theory. In the 
general framework of classical DFT crystals are viewed as "self-sustained", periodic den- 
sity oscillations of a liquid, which minimize a unique, but in general unknown free energy 
functional. Ramakrishnan and Yussouff demonstrated llj that a simple functional, which 



is Taylor-expanded about a homogeneous liquid state near coexistence semi-quantitatively 
accounts for the freezing transition in the hard sphere system. Such Taylor-expanded func- 
tionals can be devised for a wide range of two-particle potentials but they are often not very 
precise. Nevertheless they are a useful starting point for deriving more coarse-grained models 
via gradient expansions leading to phase field crystal models for materials science [12]. For 
hard-body potentials there is a constructive way to derive functionals "from scratch" (not 
relying on perturbative expansions) using essentially geometric arguments. This approach 



is known as fundamental measure theory (FMT) 13l4l5l| . With regard to the description of 



crystals, it has proved to be fruitful to consider the zero-dimensional (OD) limit of density 



distributions localized to a point and their exactly known free energy 



161 ]. By requiring 



that the density functional reproduces this OD free energy for density peaks at one, two and 



three points in space, a density functional may be constructed which exhibits sohd phase 
properties in very good agreement with simulations [17|. (In the case of density distributions 
with 5-peaks at three points, the OD free energy is reproduced only approximately.) 

In the seminal work 



17l |. the crystal density distributions were parametrized with isotropic 
Gaussians with variable width parameter and normalization (to allow for a finite vacancy 
concentration n^ac)- By minimizing the free energy with respect to the width parameter 
and the normalization, the following results were obtained: The crystal free energy per 
particle F/N agrees with simulation to within less than a percent and the Gaussian width 
is only slightly smaller than seen in simulations. However, furthermore it was found: Z(i) 



No free energy minimum for rivac > and (ii) equal free energies for 



C 



be and hep. In a study 



combining simulation and free minimization of FMT functionals [181 it was shown that (i) 
is a defect of the functional used in [17[ (the Tarazona tensor functional) and that upon 
free minimization the White Bear II tensor functional of Ref. [ij] gives thermodynamically 
consistent result^j with a small equilibrium vacancy concentration ra^ac ~ 2-10~^ for fee. The 
free energies per particle obtained by free minimization vs. constrained minimization using 
isotropic Gaussians differ by about 2 ■ 10"'^ k^T (near coexistence), which is of the order 
of magnitude one would also expect for the fcc-hcp difference AF/N. Hence one is lead to 
the suspicion that (ii) (i.e. AF/N = 0) is an artefact of the constrained minimization. This 
issue will be addressed here. 

Apart from the issue of fee vs. hep in hard spheres, FMT is also suited to investigate the 
metastable bee crystal (which in FMT is simply stabilized by the periodic boundary con- 
ditions). A previous FMT study [19|] using constrained minimization found two metastable 
bcc branches as well as a peculiar behavior of the lattice site density peaks when the density 
is increased. We will investigate this finding further by fully minimizing the FMT functional 
and will relate our results to the Stillinger series. 

The article will be structured as follows: We recapitulate basic FMT as used here 
(Sec. \llM and Stillinger's expansion in correlated, contiguous particles (Sec. IIIBl) . Results 
from both approaches are presented in Sec. IIIII and Sec. IIVI summarizes and concludes our 
work. 



^ This is discussed in Ref. 18|, Sec. III. A. under the heading "ju consistency". 



II. THEORY 

A. Fundamental measure theory 

1. Definition of functionals 

In the framework of density functional theory, the grand canonical free energy is a func- 
tional of the one-body density profile p(r) 

m = r^lp] + ^-[p] - J rfr(/i - V^^\r))p{r) . (1) 

where J-""^ and J-"*^^ denote the ideal and excess free energy functionals of the fluid, /i denotes 
the chemical potential and the external potential is represented by V^^^. The exact form of 
the ideal part of the free energy is given by 

l3J^''^[p\ = f d^r(3f'^{r) = f d^rp{r)i\n[A^p{r)] - 1) . (2) 

Here, A is the de-Broglie wavelength and /3 = l/lk^T). 

Fundamental measure theory (FMT) currently is the most precise functional for the excess 
free energy part for the hard sphere fluid. The corresponding excess free energy is given by 

:F^^ = I drr{{n[pir)]})) , (3) 

77.1722 - 111 ■ n2 



/3r{{n[p{r)]})) = no ln(l - n,) + ^,{n,) 



1-713 



3 (-772 na ■ n2 + ^2,^77*^ 7^2 j + 71277* -71*^ - 71* .7Z*fc7z|J 
+^2 ^3 —^ ^2 • ^ 

167r(l — Tis)^ 



Here, f^^ is the excess free energy density which is a (local) function of a set of weighted 
densities {n(r)} = {770,711,712,773, rii, n2, 777'} with four scalar, two vector and one tensorial 
weighted densities. These are related to the density profile p(r) by the convolutions ^^(r) = 



J dv' p{t') w"(r — r'). The weight functions are given by {R = a/2 is the hard sphere radius): 



By choosing 



w^{r) 


= e{R-r), 


2/ \ 

w (r) 


= 5{R-r), 


w\r) 


= w\v)/{AtiR), 


w\r) 


= ^^(r)/(4vri?2), 


w2(r) 


= Y/r5{R-r) , 


w^(r) 


= wV(47ri?) , 


wlj 


= riTj/r'^ 6{R-r) . 



(5) 



(/?i = 1 and v^2 = 1 



(6) 



we obtain Tarazona's tensor functional [l7| based on the original Rosenfeld functional 13 |. 
The choice 



^1 

^2 



1, 
1 - 



-2n, 



3nl 



2(1 - ri3)2ln(l - ns) 



3nl 



(7) 



corresponds to the tensor version of the White Bear I functional [20( . Finally, with 

2n3 -nl + 2(1 - n3)ln(l - ns) 



^1 



^2 



3nl 



2n3 - 3n| + 2n| + 2(1 - n3)2ln(l - Us) 



(8) 



the tensor version of the white Bear II functional is recovered [ij]. This functional is most 
consistent with respect to restrictions imposed by morphological thermodynamics 21 1. 

In density functional theory, the crystal is viewed as a self-sustained inhomogeneous 
fluid. Therefore, beside bulk and inhomogeneous fluids, it is possible to study properties 
of the hard-sphere crystal within the framework of FMT. Using the variational principle, 
the equilibrium density profile Pcq(i") is determined via minimizing the grand canonical free 
energy functional which leads to the Euler-Lagrange equation: 

Peq(r) _ 5^^^[p(r)] 



/3-Mn 



Po 



5p{r) 



^ 



V 



(9) 



For the equilibrium crystal, V^^^i^r) = and Pcq(r) is lattice-periodic, and po, the homo- 
geneous density (bulk density), is fixed by the excess chemical potential /i^^. Being com- 
putationally simpler than a free minimization of the density profile, crystal density profiles 
are often obtained by a constrained minimization of a model profile with only a few free 
parameters such as e.g. a Gaussian profile 

Pcr(r)= Yl (l-^vac)fH exp (-a{r - ri)A . (10) 

lattice sites i \ / \ / 

Here, the free parameters are the Gaussian peak width a and the vacancy concentration 

2. Choice of unit cells for the numerical solution of Euler- Lagrange equation 

Face centered cubic (fee) and hexagonal close-packed (hep) are two regular lattices with 
the highest possible hard-sphere packing fraction (77 ^ 0.74). The body centred cubic (bcc) 
structure can attain only packing fractions up to r/ ^ 0.68. The fee and hep structures differ 
in how sheets of hexagonally packed hard spheres are stacked upon one another. Relative 
to a reference layer A (see Fig. [1]), two other layer types B and C are possible which are 
laterally shifted with respect to A. In the fee structure the stacking of the hexagonally- 
packed planes corresponds to the crystallographic [HI] direction and every third layer is 
the same (ABCABCA) whereas in the hep lattice ([001] direction), the sequence of A and B 
repeats (AB ABABA) (Fig. [1]). If the binding energy (or free energy) were dependent only 
on the number of nearest-neighbor bonds per atom (bonds have no direction), there would 
be no energetic difference between the fee and hep structures. 

The most convenient unit cell for fee is the cubic unit cell with 8 particles at the corners 
and 6 face-centred particles (this cell, however, lies oblique in the ABCABCA packing dis- 
cussed above). For hep it is the unit cell with hexagonally packed hard spheres on the basal 
plane. In order to avoid any numerical errors in the comparison between fee and hep, we 
define two extended unit cells of the same size with hexagonally packed spheres as the base 
plane (see Fig. [1]). Discretizing the extended unit cells by the same number of equal-distant 
grid-points ensures that the lattice points in layer A are on grid points and for layers B and 
C the lattice points are equally "off-grid" since there is a mirror reflection symmetry with 
respect to the x-axis between B and C. In view of the narrow density peaks centered around 



each lattice point, this choice ehininates numerical differences between fee and hep free en- 
ergies to a large extent. In Fig. [H a is the nearest neighbor distance, and in the close-packed 
case a = a. The fee cubic symmetry requires c = a/8/3 a which entails that the distance 
between nearest neighbors within a base plane is the same as between neighboring planes. 
For hep, the hexagonal symmetry group does not enforce this constraint, for a discussion of 
the implications thereof see Sec. IIIICI below. 





FCC 



HCP 



FIG. 1. (color online) Extended unit cells for the fee and hep crystal structures. The fee layers 

cycle among the three identical, but laterally shifted layers, the blue A layer, the red B layer and 

the green C layer. For hep, the A and B layers alternate. Positions of the lattice points of the first 

layers from the bottom are: 

layer A: (0,0,0), (a, 0,0), (0, \/3a,0), (a,V^a,0), (^a, ^a,0), layer B : (0, ^o, ^c), 

{a,^a,^c), (ifl, ^o, ic), layer C : {^a,^a,c), (0, ^a,c), {a,-^a,c). 

a is nearest neighbor distance in the basal plane and c/2 is the distance between two neighbouring 

layers. 



3. Free minimization 

We determine the equilibrium crystal profile pcq{r; Pq, n^ac) by a full minimization in 
three-dimensional real space. The density p(r) is discretized over a cuboid volume with 
edge lengths L^ = Ly = L^ = a for bee (using the cubic unit cell) and L^ = a, Ly = y/Sa 
and Lz = 3c for both fee and hep (using the extended unit cells of Fig. [1]) with periodic 
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boundary conditions. We perform a double step minimization of the free energy. In the 
first step, the bulk density po and the vacancy concentration rivac are fixed and the Euler- 
Lagrange equation (^ is solved iteratively with a start profile given by the Gaussian profile 
f ITU]) with optimal width. The excess chemical potential fi'^^ in Eq. (^ is treated as a Lagrange 
multiplier to ensure the constraint of fixed rivac- In the next step, this procedure is repeated 
for different rzvac (still keeping po fixed), and the equilibrium density profile is determined 
by minimizing the free energy per particle with respect to the the vacancy concentration. 



nvac- For a more detailed discussion of this procedure see Ref. [18 1. 



In the program, the density profile p and 11 weighted densities (two scalar densities 
77-3, 77-2, three vector densities, U2,i, i = x,y,z, and six tensor densities, u^j) need to be 
discretized on a three dimensional grid covering the cuboid boxes. Usually we chose grids 
for the bcc unit cell with 64 x 64 x 64 points in the x, y and z directions, respectively, 
and 128 x 128 x 384 points for the fee and hep extended unit cells. Convolutions in real 
space are multiplications in Fourier space. The necessary convolutions are computed using 
Fast Fourier Transformations. We use the FFTW 3.3 library for parallelized Fast Fourier 
Transforms. The other parts of the code are parallelized through OpenMP. 

There are many sophisticated algorithms for minimizing a function and likewise many 
techniques to increase the speed and efficiency of the process. To have a more efficient 
algorithm, the iteration of Eq. (Q was done using a combination of Picard steps and DIIS 



steps (Discrete Inversion in Iterative Subspace) 22|. In order to prevent the procedure from 
diverging during the Picard iterations, in each step we mix the new density with the old 
one, 

Pnew = (1 - a)Pold + aPnew • (H) 

Here, a is a mixing parameter and it is usually a small number. For the case of bcc, a can be 
adapted in the course of the iterations in the range of a = 10"^ . . . 10^^. For fee and hep, a 
constant value for a stabilizes the iterations, with values a = 10~^ . . . lO^"'. A typical FMT 
run consisted of an initial Picard sequence with about 30 steps. Then we alternated between 
Picard sequence of 7 steps and a DIIS step (which needs another nous Picard initialization 



steps), see also Ref. [23 1 
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B. Stillinger's expansion in correlated, contiguous particles 



1. General outline 



Consider the canonical partition function for TV hard spheres: 



N 



Q{N,V,T) = j;^Jdr,...JdT^ n '^^^^■)' (12) 



i,j {i<j) 



0(zj) = ^,'- . (13) 



1 {Vij > a) 



Here, Vij = |rj — r^l is the center distance between particles i and j. We consider a reference 
lattice of our choice (fee, hep or bcc) with M > N lattice sites at positions Sj spanning the 
volume V. We associate each particle i with a lattice site at site Sj and that association 
divides the 3A^ dimensional configuration space into nonoverlapping regions Qi^p. The precise 



form of this association is discussed in Ref. 



5| , but one may think of it loosely in terms of 



each particle i belonging to the Voronoi cell around site Sj of the lattice. For a chosen subset 
of A^ lattice sites {sj} and associated cells, the index p runs over the A^! permutations of the 
particles among these cells and this leads to an identical division of the configuration space, 
Qi^pj^ = ^i,p2- The index / runs over the different associations of A^ particles with M > N 
lattice sites and becomes important in the case of finite vacancy concentration. Thus we 
obtain for the partition function: 



Q(N,V,T) = -^J2 [■■■ [ dv,...dr^ll<l,itj). (14) 



For zero vacancy concentration, this decomposition is akin to the SOC method (as e.g. 



Is 



discussed in Ref. [9[) where each particle is confined to its Wigner-Seitz cell. Following 
Ref. p], one may write Q in terms of configuration integrals Z^, Z^- , ...which describe 
the correlated motion of one, two, . . .particles in a background matrix of A^ — 1, N — 2, 
. . . particles fixed at their associated lattice sites. These configuration integrals are defined 

10 



as 

N 



Zl= f dr,ll(f){tj) with (15) 



rj = s, (j ^ i 

N 



Z\. = I dvidrj Yl 0(^^)0(jfc) with (16) 






The integration domains must fulfill u\, uj\p ■ ■ ■ G ^2/ 1, and they depend on the indices of the 
free particles i, j and also in the index / determining at which lattice sites the other particles 
are fixed. The partition function is now expressed as the product 

-1 N N yl N yl yl yl 7/ 

^^ I i i<j ^i^i i<j<k ^ij^ik^jk 

^ N N 

-A^n^' n^; u^k (18) 

i i<j i<j<k 



The Y's can also be expressed by the recursive relation 

yLn = r, ^Hi ' (19) 



where {^i . . . im} is any proper subset of {1 . . . n}. (For example, when omitting indices we 
have F2 = Z^/{Y{f2) and F3 = Z:,/{Y{f2Y:,Y,2Yi^Y2:,).) 

2. Expansion up to n = 2 for hep, fee and bee hard spheres 

In the following, we restrict calculations to the case N = M (number of particles equal 



to number of lattice sites), i.e. consider a vacancy-free crystal. From simulations [2^ and 
FMT [l8| we can estimate that the effect of vacancies on the free energy of the crystal is 
small: for fee hard spheres we have rivac ~ 10~^ (simulations) and rivac ~ 10~^ (FMT) in 
equilibrium at coexistence, the corresponding free energy shift compared to n^^c ~^ can be 
estimated from FMT, AF/N ~ 10"^ ksT. 

Truncated after the first term, the Stillinger series is 

^1 = J^(V,f , (20) 

11 



where Z\ has been reduced to V\ , the free volume for one particle in a cage of fixed neighbors 
at their lattice sites. Consequently the free energy is 

/3Fi = -iVln^. (21) 



For fee and hep, V\ is equal and has been calculated analytically in Ref. [25|, we quote 
this result in App. |Al For bcc, we did not find a literature result and therefore give the 
calculation and result also in App. |Al 

The second term in the Stillinger series for Q gives only a contribution different from 1 
if the two fixed particles are neighbors. Thus the truncated Stillinger series is 

* = X^(^'.)"n(if5) (22) 

Here, V2,fc is the correlated free volume of the two neighboring particles (with dimension 
(length)^) which may depend on the type of neighbor configuration (index k\ The power 
g^N reflects the freedom to choose the first of the two particles to be any of the A^ particles 
in the system and g^ is the multiplicity of the neighbor configuration. It is half the number 
of neighbors of type k for a given fixed particle. The associated free energy is 

/5F2 = /5Fi - AT ^ (7fc In {^^ ■ (23) 

For our considered lattice cases the neighbor types and multiplicities are given in Tab. [11 The 
cubic lattices fee and bcc have only one neighbor type whereas for hep there is a difference 
whether the neighbor is within the same close-packed plane or in an adjacent close-packed 
plane. See also Ref. |7| for the multiplicities corresponding to the third term in the series 
(fee and hep). 

We calculate the two-particle volumes V^2,fc for different densities by a simple Monte- 
Carlo computation. For that we specify a suitably large cuboid volume Vc for each of the 
two free particles from which n sets of random positions (for each of the two particles) are 
drawn. For each set of random positions overlap is checked with the other particle and 
the fixed neighboring particle, leading to a total of v! sets of random positions with no 
overlap. Then V2,fc = (n'/n)\/?. The statistical error Ay2,fc/V2,fc needs to be below 10~^ for 
a reliable assessment of the free energy difference between fee and hep, and this is achieved 
with 1000 subsets, each containing n = 10^ sets of random positions. In the limit po — ^ Pep 
(pcp = Y^/o"^ is the close-packing density) agreement was found with the analytical results 
of Ref. [TJ, but we had to approach p^p very closely to establish that. 

12 



lattice neighbor type k g^ 

fee all neighbors 1 6 

hep within close-packed plane 1 3 

in adjacent close-packed planes 2 3 

bcc all neighbors 1 4 

TABLE I. Neighbor configurations with multiplicities for the different lattices. 



III. RESULTS 



A. Stillinger series 



For fee and hep, the Stillinger series truneated at n = 2 gives very good results for the 
free energy per particle F/N (see Fig. |2l to obtain numbers, we put A = a). We have 



compared to very precise simulation data obtained in Refs. [18|, |26| which have an error of 
about 0.002 k^T. The Stillinger series {n = 2) results for F/N deviate from these ranging 
from 0.01 ksT (at poa^ = 1.0) to 0.03 ksT (at poa^ = 1.15), this is less than 0.5% relative 



deviation. This is about the same accuracy we obtain with FMT (see also Ref. [18[). Note, 
however, that a deviation of the order of 0.01 k-^T is about 10 times higher than the fcc-hcp 
free energy difference obtained from simulations, as discussed before. 

For bcc, the situation is very much different. Since the bcc structure for hard sphere is 
unstable against shear, the crystal can be stabilized in simulations only by constraints such 
as in the SOC method. We would expect from the previous derivation that the Stillinger 
expansion is a reasonable series expansion for the free energy of the SOC method. However, 
as Fig. |2] demonstrates, the first two terms are quite far away from the SOC data and also 
from the FMT results for the branch with lowest free energy, pointing to the importance 
of higher correlations. (Ultimately, the shear instability is a collective many-body effect, 
so perhaps the importance of many-particle correlations also in the constrained crystal is 
not too surprising.) See, however, the next subsection for a more detailed discussion on bcc 
solutions within FMT, especially with regard to a solution branch with higher free energy 
which appears to be linked to the bcc Stillinger solution. 

13 
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FIG. 2. (color online) Crystal free energies J3F/N for fee and bcc from the Stillinger series in 
comparison to simulation data and FMT results (bcc). For fee, simulation data are taken from 



Refs. 
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26l |. and for bcc, simulation data are obtained using the single-occupancy cell method 



(SOC) [27J. The FMT data are this work, see Sec. 111111 



Finally, for fcc/hcp the inclusion of the correlated neighbor term increases the free energy, 
whereas for bcc it leads to a decrease. 



B. bcc FMT results 

As already discussed, a bcc crystal solution can only be stabilized by constraints. In 
FMT, these are the periodic boundary condition on the cubic unit cell. Within the Gaussian 
parametrization (see Eq. (fTOj) ). bcc solutions in FMT (Rosenfeld, Tensor and White Bear 
Tensor, see Sec. IIIAp have been investigated by Lutsko 19| (with the additional constraint 
'^vac = 0, such that in the free energy minimization, the width parameter a is the only 
variable which is varied at a given bulk density po)- For small bulk densities {poc^ ^ 
1.16), Lutsko found a single free energy minimum with a rather small width parameter 
a ~ 30 . . .40, indicating a broad Gaussian peak. Interestingly, a(po) exhibits a maximum 



at poc^ ~ 1-13 and then decreases again upon increasing the density (i.e. the density peaks 
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FIG. 3. (color online) (a) Difference in free energy per particle between the fully minimized and 
the Gaussian solution for the first branch of the bcc solutions as a function of bulk density. Inset: 
Equilibrium vacancy concentration as a function of bulk density for the same first branch, (b) Free 
energy per particle as a function of bulk density for the bcc solution of the second branch: Full 
minimization (symbols, n^ac = 6 x 10~^ fixed) and Gaussian approximation (full black line). For 
comparison the Stillinger result (n = 2) is given (dashed line) as well as the Gaussian approximation 
for the first branch (dot-dashed line). 

become wider upon compressing the crystal!). Moreover, at bulk densities Pqct^ > 1.16 a 
second free energy minimum was visible (with higer free energy). In this second branch, the 
width parameter increased (the peak width decreased) with increasing density as one would 
naively expect. 

We investigate these findings further using full minimization. For the first branch with 
lowest free energy, we confirm that there is a minimal width of the peaks at poa^ ~ 1.13. 
Full minimization reveals a rather strong deviation from the simple Gaussian form in the 
density peaks: The difference in free energy per particle F/N between Gaussian and full 
minimization is about 0.1 k-^T (see Fig. |3] (a)) and thus about 2 orders of magnitude higher 



than in the case of fee [18j . Curiously, this free energy difference increases with increasing 
density beyond poa^ ~ 1.07. Secondly, the equiUbrium vacancy concentration n^ac is of 
the order of 10"^ and thus several orders of magnitude higher than found in fee. n^acipo) 
has a minimum at p^a^ f^ 1.10 and then increases again, adding to the peculiarities of this 
solution branch. We note that in an FMT study of parallel hard squares and cubes similar 
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peculiarities have been found [28 1. 



The second branch found by Lutsko is not an artefact of the constrained Gauss minimiza- 
tion. By a careful iteration procedure, we found corresponding fully minimized solutions 
whose free energy per particle is very close to the values from the Gaussian approximation 
(thus very much like the fee solutions and very much unlike the solutions from the first 
branch), see Fig. 12] (b). For increasing densities, we see a convergence of F/N to the re- 
sults of the Stillinger series {n = 2). Thus the second branch of the bcc solutions has the 
same character as the fee solution when compared with the Stillinger approach: only a few 
correlated particles are sufficient to obtain the free energy. 

One could argue that the discussion of these bcc solutions is futile and void of physical 
significance in view of their overall instability. However, the quality of the FMT functionals 
and their success in describing the fee phase leads us to think that these solutions are perhaps 
not to be discarded altogether. Since around coexistence {pocr^ ~ 1.04) the difference in F/N 
to the fee crystal is about 0.3 k-oT and thus very high, it is reasonable that bcc crystallites 
have not been observed in the nucleation process of a hard sphere crystal. Nevertheless, the 
bcc solutions are perhaps a useful reference point for discussing the crossover from fee to 
bcc as the most stable crystal structure for other potentials such as of (a/r)" type. These 
could be treated by suitable perturbation ansatz in the free energy functional. Also, it could 
be interesting to investigate further the dispersion relation of phonons for the solutions of 
the first branch and thus shed further light on the shear instability. 



C. fcc/hcp: Free energy differences and density anisotropies 



As discussed in the Introduction, FMT gives the same free energy per particle F/N for 



fee and hep when the Gaussian approximation is employed [17|, ll9| . Free minimization lifts 
this degeneracy in the free energy. In order to understand this result qualitatively, it is useful 
to consider the symmetries in the unit cell of fcc/hcp and the constraints these symmetries 
place upon the lattice-site density profiles. For fee, this is best discussed by considering the 
cubic unit cell in Fig. H] (al). The non-radial contributions to the density profile around the 
lattice point in the origin can be expanded in a Taylor series in x,y, z where the terms in 
this series must respect the 48 point symmetry operations in the cubic unit cell (belonging 
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FIG. 4. (color online) Unit cells and density anisotropics for fee and hep. (al) and (a2) show the 
most eonvenient unit cells (cubic for fee and hexagonal for hep) for the mathematical discussion 
of the density anisotropics (see Eqs. (j24p and (j25p ). (a3) and (a4) show the unit cells used in the 
numerical computations. The hexagonally packed planes (marked in different colors) lie oblique 
in the cubic unit cell (al). (b) fee and hep density distributions around the lattice site at the 
origin in different directions. Here, we used the bulk density poa"^ = 1.04 and fixed the vacancy 
concentration to nvar = 10~^. 



to point group —3— in Hermann-Mauguin notation) 181 ]: 



Piccix, y, z) = pri,d{r) (1 + i^4(a;^ + 2/^ + 2^) + 



)• 



(24) 



Here, Pradl'") is an averaged, radial profile which is more or less of Gaussian shape. The 
leading anisotropic term is of polynomial order 4 with expansion coefficient K4. One can 
also understand this result by resorting to an expansion in the subset of spherical harmonics 
which respect the cubic point symmetry, this leads to an expansion in the so-called Kubic 
Harmonics 29||. - For hep, we consider the unit cell in Fig. H] (a2). The corresponding 



Taylor expansion for the non-radial contributions to the density profile around the lattice 
point in the origin has to respect only the 24 point symmetry operations appropriate for the 



hexagonal group 



m m m 



. According to Ref. |30|, this leads to 



Phcp{x, y, z) = prad(r) (1 + K'^z + K'^yiSx -y ) + ...) , 



(25) 



where polynomial terms up to order 3 have been taken into account (with expansion coeffi- 
cients K^). The corresponding construction using spherical harmonics leads to the so-called 
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Hexagonal Harmonics. We observe that there is a quahtative difference in the shape of the 
density profile between hep and fee according to these expansions: 

(z): To leading order in anisotropy for hep, the density peak p{r) should look different in z- 
direction (perpendicular to the hexagonally packed planes) than in directions in the x- 
y plane. To phrase it differently: one would expect different width parameters a^, ax,y 
for a Gaussian density peak of the form phcp{x, y, z) ex: exp(— ^^^^^(x^ + y"^) — azZ^). We 
did not observe this in our numerical solutions but we will return to this point below. 

{ii): To next-to-leading order in the anisotropy for hep, we expect a different behavior 
when comparing p(0, y, 0) with p(0, —y, 0) due to the antisymmetric term oc K'^ in 
Eq. ( 12 5p . Such a symmetry breaking is not present in the fee peak. To demonstrate 
this difference, we cornpare p{0,±y,0), p(x,0,0), and p{z, 0,0) between fee and hep, 
see Fig. H] (b) and (c)o Indeed we observe that the symmetry is broken for the hep 
profile, in accordance with the anisotropy expansion, and we conclude that the fcc/hcp 
free energy difference in FMT results from this symmetry breaking. 

Our results for the fcc/hcp free energy difference per particle are given in Fig. [5t^a). In 
FMT (White Bear H-Tensor), the difference fiAF/N is larger than zero, implying that hep 
has lower free energy. Furthermore, there is only a moderate drop of /3AF/N with the 
bulk density po- At coexistence (poc'^ = 1-04), we have computed /3AF/N also for other 
FMT functionals (Tarazona-Tensor, White Bear-Tensor) and found no change in sign but a 
variation in magnitude by 50% or 5 ■ 10"^. In view of the variation of I3F/N for fee between 
the functionals (about 4 ■ 10"^, i.e. a factor of 80 larger), the functionals are very consistent 
with each other with respect to the stability of hep. The results from the Stillinger series 
{n = 2) for PAF/N are approximately constant (~ 1 ■ 10^'^) with increasing density and 
coincide with the analytical value at close packing obtained in Ref. |7|. It is remarkable that 
also the FMT results seem to converge to this value. - For comparison, in Fig. |5t^a) we have 
also included the analytical value from the Stillinger series (n = 5) [Sj and the simulation 
^ Note that in our numerical computations we used the unit cells depicted in Fig. |4] (a3) (fee), and in 
Fig. |4] (a4) (hep). Thus, the fee cubic unit cell and the unit cell in Fig. [4] (a3) are related by a three- 
dimensional rotation. Likewise, the anisotropy expansion for the extended unit cell must be obtained 
from the corresponding expression ()24p for the cubic unit cell by applying this rotation. However, since 
the density anisotropy is oc y^ (a; = 0, z = 0) in Eq. ()24|) . the corresponding density anisotropy must also 

be ex y'^ {x' = 0, z' = 0) in the rotated unit cell (primes denote the coordinates in the extended unit cell 
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inFig.|l(a3)). 



value of Ref. [10|. Although FMT does not agree with the sign of PAF/N obtained in the 
simulation, it is gratifying to note that according to these results FMT is correct on the 
level of two correlated particles in the Stillinger picture. 

Finally, we return to the observation that in the hep density anisotropy the leading term 
oc z^ (see Eq. (l25l) ) was missing in our numerical solutions. This is related to our choice 
of the distance between the hexagonally packed layers (c/2 = cq/2 = ^/2/3a where a is 
the nearest neighbor distance, see Fig. [1]). With this choice the distance between nearest 
neighbors is the same for two sites within the same hexagonally packed planes and two sites 
in two adjacent planes. However, the hep symmetry group does not require this, and one 
is free to choose another distance between the planes. With a different choice, also the 
nearest neighbor distance is different for sites in two different planes and also the width of 
the lattice site density profiles will be different in the direction normal to the hexagonally 
packed planes. We have investigated whether also the free energy minimum for hep shifts to 
a value different from cq. In order to keep the bulk density constant we defined a stretching 
parameter, 7 = c/cq, which describes the distortion of the crystal in 2;-direction. In order 
to keep the bulk density constant, we rescaled the nearest neighbor distance in the planes 
as follows: a' = 0/^/7. Full minimization was done for a range of 7 values. The result for 7 
which minimizes F/N is shown in Fig. |5] and it is seen that the equilibrium distortion is quite 
small, below 10~^. The corresponding free energy shift per particle compared to the solution 



with c = Co is about 10"^ A;bT. These results are actually similar to the ones in Ref. 31| : 
There, a similar lattice distortion was calculated for the zero-temperature Lennard- Jones 
hep crystal by lattice sums. 

IV. SUMMARY AND CONCLUSIONS 

In this work we have performed a study of bcc, fee and hep hard sphere crystals using 
unrestricted minimization in density functional theory (DFT) of Fundamental Measure type 
(FMT) which is currently the most accurate approach. We have complemented these inves- 
tigations with an approach which is based on the expanding the crystal partition function 
in terms of number n of free particles while the remaining particles are frozen at their ideal 
lattice positions (Stillinger series). 

For the metastable bcc crystal, we have found two solutions for bcc crystals whose free 
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FIG. 5. (color online) (a) Free energy difference between fee and hep vs. bulk density. The black 
symbol shows the simulation value from Ref. [lOt]. Rest of the symbols show the data obtained 
from FMT and the Stillinger series (n = 2) and dashed lines show the asymptotic behavior of the 
free energy difference near close packing for the Stillinger series (different n) [8[. (b) Distortion 
parameter 7 = c/cq which mimimizes the hep free energy vs. bulk density. In all FMT calculations 
we put rivac = 10~^. 



energies are well above the free energies of fcc/hcp (see Fig. [2]and[3](b)). The first solution 
(with a rather large density peak width at lattice sites) is characterized by a rather large 
equilibrium vacancy concentration (~ 0.01) and its free energy can not be described by the 
Stillinger approach. The shear instability of bcc is presumably related to this first solution. 
The second solution (characterized by a small peak width and small equilibrium vacancy 
concentrations) agrees well with the solution from the Stillinger approach [n = 2) with 
respect to its free energy. 

The free energy degeneracy between fee and hep, found in previous approaches using 
constrained, rotationally-symmetric density peaks around lattice sites, is broken upon full 
minimization. The density asymmetries are qualitatively different for fee and hep and agree 
with expansions in respective lattice harmonics (see Fig. HJ. We found that in FMT the free 
energy per particle is lower for hep than the one for fee by about 10"'^ k^T. This agrees 
remarkably well with the Stillinger solution for n = 2 (see Fig. \5^. Simulations, however, 
indicate that fee has a lower free energy than hep by about the same figure. Previous 
investigations of the Stillinger approach in the high-density limit (near close packing) have 
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shown that hep is more stable than fee for n = 2 . . . 4 and the situation reverses for n = 5. 
Thus, the stabihty of fee seems to be a subtle effect involving the correlated motion of at 
least 5 particles which currently can not be captured by the FMT functionals. 

Appendix A: One particle volumes for the fcc/hcp and bcc hard sphere crystal 
1. fee and hep 



The one-particle free volume is equal for fee and hep and has been given in Ref . 25| . We 
introduce the nearest neighbor distance d = 2^/^Pq . The hard sphere diameter is a and 
the formula is valid for densities poa^ E [1/2, -\/2] : 






2v2(c^ — 6ccr^) I arcsin — h arcsinm ) + (Al) 

\ Q J 



la'^ { 2 arcsin u -\ arcsin w — arcsin t 



with 



c 



d/V2 



s = V3a2 - 2c2 



q = v2o^— (? • 

m = {c — 2s)/(3g) , 

t= {a^ + ca- c^)/{qa) , 

u = [{2a + c){a+ [2c - s]/3) - {a + c)^]/[q{cT + [2c - s]/3)] 

w = (cr^ — ca — (?)/{qa) , 



The shape of the free volumes is sketched in Fig. [61 



bcc 



In case of bcc the free volume is given by an octahedral-like body (see Fig. E]) centered in 
the cubic unit cell. The faces are parts of the surfaces of the exclusion spheres (of radius a) 
around the corners of the cubic unit cell. Let a = (2/po)^^'^ be the side length of the cubic 
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FIG. 6. (color online) Shape of one-particle free volumes for fee, hep and bee (from left to right) 
at a crystal density oi p^a^ = 1. 

unit cell. The free volume is then given by 

(A2) 




+ a 



ir'-r'){ 



2c vr 

arctan 

a 4 



a"' 1 n a^ c 

— c H — a arctan arctan — 

4 3V 4ac a 



(A3) 



c= V^' - o?/2 . 
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